package com.alibaba.alink.common.linalg.jama;

/**
 * Cholesky Decomposition.
 * <P>
 * For a symmetric, positive definite matrix A, the Cholesky decomposition
 * is an lower triangular matrix L so that A = L*L'.
 * <P>
 * If the matrix is not symmetric or positive definite, the constructor
 * returns a partial decomposition and sets an internal flag that may
 * be queried by the isSPD() method.
 */

class CholeskyDecomposition implements java.io.Serializable {
	private static final long serialVersionUID = -693958715904680364L;

/* ------------------------
   Class variables
 * ------------------------ */

	/**
	 * Array for internal storage of decomposition.
	 *
	 * @serial internal array storage.
	 */
	private double[][] L;

	/**
	 * Row and column dimension (square matrix).
	 *
	 * @serial matrix dimension.
	 */
	private int n;

	/**
	 * Symmetric and positive definite flag.
	 *
	 * @serial is symmetric and positive definite flag.
	 */
	private boolean isspd;

/* ------------------------
   Constructor
 * ------------------------ */

	/**
	 * Cholesky algorithm for symmetric and positive definite matrix.
	 *
	 * @param Arg Square, symmetric matrix.
	 * @return Structure to access L and isspd flag.
	 */

	public CholeskyDecomposition(JaMatrix Arg) {

		// Initialize.
		double[][] A = Arg.getArray();
		n = Arg.getRowDimension();
		L = new double[n][n];
		isspd = (Arg.getColumnDimension() == n);
		// Main loop.
		for (int j = 0; j < n; j++) {
			double[] Lrowj = L[j];
			double d = 0.0;
			for (int k = 0; k < j; k++) {
				double[] Lrowk = L[k];
				double s = 0.0;
				for (int i = 0; i < k; i++) {
					s += Lrowk[i] * Lrowj[i];
				}
				Lrowj[k] = s = (A[j][k] - s) / L[k][k];
				d = d + s * s;
				isspd = isspd & (A[k][j] == A[j][k]);
			}
			d = A[j][j] - d;
			isspd = isspd & (d > 0.0);
			L[j][j] = Math.sqrt(Math.max(d, 0.0));
			for (int k = j + 1; k < n; k++) {
				L[j][k] = 0.0;
			}
		}
	}

/* ------------------------
   Temporary, experimental code.
 * ------------------------ *\

   \** Right Triangular Cholesky Decomposition.
   <P>
   For a symmetric, positive definite matrix A, the Right Cholesky
   decomposition is an upper triangular matrix R so that A = R'*R.
   This constructor computes R with the Fortran inspired column oriented
   algorithm used in LINPACK and MATLAB.  In Java, we suspect a row oriented,
   lower triangular decomposition is faster.  We have temporarily included
   this constructor here until timing experiments confirm this suspicion.
   *\

   \** Array for internal storage of right triangular decomposition. **\
   private transient double[][] R;

   \** Cholesky algorithm for symmetric and positive definite matrix.
   @param  A           Square, symmetric matrix.
   @param  rightflag   Actual value ignored.
   @return             Structure to access R and isspd flag.
   *\

   public CholeskyDecomposition (Matrix Arg, int rightflag) {
      // Initialize.
      double[][] A = Arg.getArray();
      n = Arg.getColumnDimension();
      R = new double[n][n];
      isspd = (Arg.getColumnDimension() == n);
      // Main loop.
      for (int j = 0; j < n; j++) {
         double d = 0.0;
         for (int k = 0; k < j; k++) {
            double s = A[k][j];
            for (int i = 0; i < k; i++) {
               s = s - R[i][k]*R[i][j];
            }
            R[k][j] = s = s/R[k][k];
            d = d + s*s;
            isspd = isspd & (A[k][j] == A[j][k]); 
         }
         d = A[j][j] - d;
         isspd = isspd & (d > 0.0);
         R[j][j] = Math.Sqrt(Math.max(d,0.0));
         for (int k = j+1; k < n; k++) {
            R[k][j] = 0.0;
         }
      }
   }

   \** Return upper triangular factor.
   @return     R
   *\

   public Matrix getR () {
      return new Matrix(R,n,n);
   }

\* ------------------------
   End of temporary code.
 * ------------------------ */

/* ------------------------
   Public Methods
 * ------------------------ */

	/**
	 * Is the matrix symmetric and positive definite?
	 *
	 * @return true if A is symmetric and positive definite.
	 */

	public boolean isSPD() {
		return isspd;
	}

	/**
	 * Return triangular factor.
	 *
	 * @return L
	 */

	public JaMatrix getL() {
		return new JaMatrix(L, n, n);
	}

	/**
	 * Solve A*X = B
	 *
	 * @param B A Matrix with as many rows as A and any number of columns.
	 * @return X so that L*L'*X = B
	 * @throws IllegalArgumentException Matrix row dimensions must agree.
	 * @throws RuntimeException Matrix is not symmetric positive definite.
	 */

	public JaMatrix solve(JaMatrix B) {
		if (B.getRowDimension() != n) {
			throw new IllegalArgumentException("Matrix row dimensions must agree.");
		}
		if (!isspd) {
			throw new RuntimeException("Matrix is not symmetric positive definite.");
		}

		// Copy right hand side.
		double[][] X = B.getArrayCopy();
		int nx = B.getColumnDimension();

		// Solve L*Y = B;
		for (int k = 0; k < n; k++) {
			for (int j = 0; j < nx; j++) {
				for (int i = 0; i < k; i++) {
					X[k][j] -= X[i][j] * L[k][i];
				}
				X[k][j] /= L[k][k];
			}
		}

		// Solve L'*X = Y;
		for (int k = n - 1; k >= 0; k--) {
			for (int j = 0; j < nx; j++) {
				for (int i = k + 1; i < n; i++) {
					X[k][j] -= X[i][j] * L[i][k];
				}
				X[k][j] /= L[k][k];
			}
		}

		return new JaMatrix(X, n, nx);
	}
}

